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ABSTRACT 

It  has  been  know:  for  quite  a while  that  if  a function  u(x,y)  is 
harmonic  in  a region  with  reentrant  corners,  there  are  almost  certainly  infinite 
discontinuities  of  the  first  derivative  of  u in  the  neighborhood  of  the 
reentrant  corner  (or  corners) . Simple  examples  are  for  an  L-shaped  region  or 
T-shaped  region.  Some  instances  of  these  have  been  treated  by  conformally 
mapping  the  region  into  the  interior  of  a rectangle.  Attempts  to  solve  the 
problem  as  first  posed  by  a finite  difference  scheme  or  a finite  element  scheme 
will  usually  give  poor  approximations  near  any  reentrant  corner  because  the 
finite  differences  or  finite  elements  have  large  truncation  errors  when  a first 
derivative  is  infinite.  When  conformal  mapping  is  tried,  the  conformal  maps 
are  usually  only  approximate,  and  similar  errors  arise,  for  more  or  less 
similar  reasons. 

In  view  of  recent  work  giving  convergent  expansions  for  u in  the 
neighborhood  of  reentrant  corners  (see  "Calculation  of  Potential  in  a Sector, 
Part  I,”  by  J.  Barkley  Rosser,  MRC  TSR  #1535)  one  can  now  give  accurate  solu- 
tions for  such  problems.  Some  experiments  with  such  regions  are  reported. 

AMS  (MOS)  Subject  Classification:  65N99 

Key  Words:  Partial  differential  equations,  Boundary  value  problems. 
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SIGNIFICANCE  AND  EXPLANATION 

There  are  numerous  situations  in  which  the  behavior  is  given  by  a 
function  which  is  harmonic  in  a region;  temperature  in  the  region,  fluid  flow 
through  the  region,  diffusion  through  a permeable  medium  in  the  region,  to 
name  a few.  If  there  are  no  obstacles  in  the  region,  there  are  numerous 
techniques  which  give  the  desired  harmonic  function  efficiently.  However, 
there  are  getting  to  be  more  and  more  cases  in  which  the  medium  has  obstacles 
in  it,  commonly  obstacles  with  corners.  Typically  a corner  of  an  obstacle 
will  be  a reentrant  corner  as  far  as  the  medium  is  concerned.  That  is,  one 
must  traverse  more  than  180°  THROUGH  THE  MEDIUM  to  get  from  one  edge  of 
the  angle  to  the  other.  It  was  discovered  more  than  twenty  years  ago  that, 
at  such  a reentrant  corner,  the  harmonic  function  that  one  is  seeking  usually 
has  infinite  derivatives.  In  such  case,  the  familiar  techniques  for  finding 
the  harmonic  function  give  very  poor  approximations  near  the  corner.  The 
present  paper  gives  a new  technique  which  will  give  high  accuracy  approxima- 
tions to  the  desired  harmonic  function  without  undue  labor.  Part  I,  which 
is  presented  here,  gives  the  theoretical  background  of  the  technique.  Part  II, 
which  \.ill  follow,  will  give  numerical  examples,  to  show  how  the  technique 
works  in  practical  situations. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  author  of  this  report. 


HARMONIC  FUNCTIONS  ON  REGIONS  WITH  REENTRANT  CORNERS,  PART  I 


J.  Barkley  Rosser 


1 . Background . A function  u(x,y)  is  said  to  be  harmonic  if 


We  approximate  (1.1)  by  the  familiar  difference  formula 


The  error  involves  fourth  derivatives  of  u.  We  suppose  that  these  exist 
and  are  reasonably  well  behaved.  This  is  usually  the  case,  so  for  a long  time 
it  has  been  customary  to  seek  an  approximate  solution  for  (1.1)  by  solving  the 
set  of  linear  equations  resulting  from  using  various  values  of  (x,y)  in  (1.2) 
If  h is  small,  so  that  the  error  in  (1.2)  is  small,  the  solution  of  (1.2) 
gives  good  approximations  to  the  values  of  u at  a set  of  grid  points. 


Unfortunately,  if  h is  smaj.1,  then  one  has  a very  large  number  of 
linear  equations  to  solve,  and  the  labor  of  computation  is  very  great.  Until 
the  advent  of  the  computer,  one  compromised  by  using  a fairly  large  value  of 
h,  to  curtail  the  calculation,  but  one  had  to  be  content  with  not  a very  good 
approxi ma tion . 


One  effort  to  use  computers  to  improve  this  situation  is  embodied  in  item 
[1]  of  the  Bibliography,  by  Kantorovich,  Krylov,  and  Chernin.  If  one  has 
values  for  u prescribed  on  the  boundary  of  a rectangle,  the  tables  in  II] 
allow  one  to  get  fairly  quickly  the  solution  of  (1.2)  inside  the  rectangle. 
(This  can  now  be  done  more  quickly  by  means  of  the  Fast  Fourier  Transform,  so 
that  [1]  is  now  obsolete.) 


To  show  the  effectiveness  of  their  tables,  the  authors  of  [1]  undertook  to 
find  the  solution  of  (1.2)  inside  an  L-shaped  region  (see  Fig.  1).  They 
prescribed  values  for  u around  the  boundary,  and  used  the  Schwarz  alternating 
procedure.  Specifically,  they  first  guessed  values  along  CF.  Using  these 
with  the  boundary  conditions,  the  tables  gave  the  solution  of  (1.2)  inside  the 
rectangle  ABFG.  In  particular,  they  gave  values  along  HC.  Using  these  with 
the  boundary  conditions,  the  tables  gave  the  solution  of  (1.2)  inside  the 
rectangle  HDEG.  This  led  to  a better  guess  for  the  values  along  CF.  Using 
the  better  guess,  the  process  was  repeated.  After  a modest  number  of  repeti- 
tions, the  procedure  converged  to  give  about  six  decimal  accuracy. 
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2.  Two  lemmas.  To  get  some  comprehension  of  the  results  of  Wasow  and 
Lehman  in  12]  and  [3] , we  use  two  lemmas,  which  we  state  here  without  proof 


Let  f (a)  = f(b)  = 0.  Let  f(x)  have  almost  every 
a < x < b which  is  of  bounded  variation.  Then 
in  the  interval  a < x < b, 


with 


can  be  calculated  very 


By  "converges  very  rapidly"  is  meant  that  |u  | goes  to  zero  at  least  of 
-3  D’ 

the  order  of  m . Thus  one  can  truncate  the  series  on  the  right  of  (2.1)  after 
500  terms  and  reasonably  expect  to  get  from  six  to  eight  significant  decimal 
places  correct.  And  the  Fast  Fourier  Transform  will  enable  one  to  calculate  the 
needed  500  coefficients  very  quickly.  The  reasoning  to  establish  this  lemma  is 
given  in  pp.  6-8  of  14) . 


Lemma  2.2.  Let  u(x,y)  be  harmonic  in  a region,  part  of  the  boundary  of 
which  is  a straight  line  segment.  Let  u(x,y)  * f(s)  on  this  straight  line 
segment,  where  s is  length  along  the  segment.  Let  f (s)  and  its  first  n 
derivatives  be  continuous,  and  let  the  (n  + l)-st  derivative  be  bounded  and 
continuous  except  at  a set  of  points  of  measure  zero.  Then  each  partial  deriva 
tive  of  u(x,y)  of  order  n has  a continuous  extension  to  the  straight  line 
boundary. 


The  truth  of  the 


In  the  present  report,  we  shall  confine  our  attention  to  the  case  where  the 
prescribed  values  of  u around  the  boundary  are  quite  smooth;  say  that  the  third 
derivative  is  bounded  and  continuous  except  at  a set  of  measure  zero.  It  is 
planned  to  write  a sequel  to  15)  explaining  how  to  handle  a variety  of  irregular- 
ities along  the  boundary.  Certain  sorts  of  irregularities  that  could  occur  along 
the  boundary  can  be  "removed"  by  the  methods  given  on  pp.  221-222  of  [6) . So  it 
does  not  seem  unduly  restrictive  to  confine  our  attention  to  harmonic  functions 
u(x,y)  which  are  quite  smooth  around  the  boundary.  In  the  present  report  we  do  so 


Let  u(x,y)  be  such  a harmonic  function  in  the  L-shaped  region  of  Fig.  2 


Indeed  we  will 
we  will  have 

(2.3) 


shortly  specialize  to  prescribing  that  on  the  boundary  in  Fig. 
u(x,y)  ■ ^ £n{ (x  + l)2  + (y  - l)2)  . 
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In  Fig.  2 we  have  shown  three-quarters  of  a circle  of  radius  A and  center 
at  the  origin.  We  undertake  to  determine  the  behavior  of  u(x.y)  inside  the 
three-quarters  circle. 

Choose  u<S)(x,y)  a function  which  in  the  interior  of  the  figure  is  harmonic 
in  the  neighborhood  of  the  sides  BC  and  CD,  including  the  three-quarters 
circle  (one  can  take  A quite  small  if  need  be) , and  which  takes  the  same  values 
along  BC  and  CD  that  are  prescribed  for  u(x,y).  Instructions  for  finding 

such  a function  u<S) (x,y)  are  set  forth  in  (SI.  If  we  have  prescribed  the 
particular  boundary  conditions  (2.3)  for  u(x,y),  such  a function  is 


(2.4) 


T «(3  £n(z  + 1 - i)  + tn(z  + 1 + i) 
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+ £n(z  - 1 - i)  - £n(z  - 1 ♦ i) } , 


where  we  have  taken 

(2.5)  z = x + iy  . 

It  can  easily  be  verified  that  the  right  side  of  (2.4)  satisfies  (2.3)  along 
the  entire  x-axis  and  the  entire  y-axis.  Also  (2.4)  is  harmonic  except  at  the 
four  points  z - ± 1 ± i.  To  keep  the  three-quarters  circle  inside  the  region 
where  (2.4)  is  harmonic,  it  suffices  to  take  A < ^2.  We  choose  such  a value 
for  A,  and  proceed. 

( s ) 

Since  u(x,y)  and  u (x,y)  take  the  same  values  along  BC  and  CD,  we 
conclude  that  along  these  two  segments 

(2.6)  u(x,y)  - u(s) (x,y) 
is  zero. 

Measure  the  angle  0 as  usual,  counterclockwise  around  the  origin  from  CD. 
On  the  three-quarters  circle,  we  will  have  (2.6)  a function  of  0 only,  since 
we  have  fixed  A.  Call  this  f(0).  As  u(x,y)  and  u(s'  (x,y)  take  values 
along  BC  and  CD  that  are  differentiable  an  infinite  number  of  times,  it 
follows  from  Lemma  2.2  that  f (0)  is  infinitely  differentiable  as  0 approaches 
0+  and  (3ir/2)-.  For  0 < 0 < 3ir/2,  f (0 ) is  the  difference  of  two  harmonic 
functions,  and  hence  is  infinitely  dif ferentiable.  So  for  0 < 6 < 3n/2,  f ( 0 ) 
has  derivatives  of  all  orders. 


As  (2.6)  is  0 on  BC  and  CD,  we  have  f(0)  = f (3n/2)  = 0.  So  we  take 
a = 0 and  b = 3tt/2  in  Lemma  2.1,  and  conclude  that  f(0)  has  a rapidly  converg- 
ing Fourier  series  expansion  for  0 0 <_  3ir/2;  put  0 for  x in  (2.1)  . 

Because  of  the  rapid  convergence,  we  see  that 
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(2.7) 


O 


I 

m*l 


sin 


2*0 

3 


is  a harmonic  function  for  0 < r < A and  0 < 0 < 3*/2.  It  equals  f(0)  for 
r * A,  and  is  0 for  6*0  or  0 ■ 3x/2.  But  (2.6)  satisfies  these  same 
conditions.  As  a h trmonic  function  is  uniquely  determined  in  a region  by  its 
values  around  the  boundary,  we  must  have  (2.6)  equal  to  (2.7)  inside  and  on  the 
three-quarters  circle.  Solving  for  u(x,y),  we  must  have  u(x,y)  equal  to  the 
sum  of  (2.4)  and  (2.7)  inside  and  on  the  three-quarters  circle. 


For  the  particular  boundary  conditions  which  we  have  chosen  (see  (2.3))  we 
have  D * 0 in  (2.7).  So  if  we  fix  a value  of  0,  0 < 0 < 3*/2,  and  approach 
the  origin  along  that  ray,  (2.7)  will  have  an  infinite  first  derivative.  As  (2.4) 
is  harmonic  inside  the  entire  circle  of  radius  A and  center  at  the  origin  (we 
took  A < Jt.)  , it  has  well  behaved  derivatives  of  all  orders  at  C.  So  u(x,y) 
must  have  an  infinite  derivative  as  r approaches  zero. 

Of  course,  if  it  had  turned  out  that  D * D * D.  ■ D&  ■ 0,  then  u(x,y) 
would  have  had  well  behaved  fourth  derivatives.  But  it  did  not  turn  out  that 
way.  In  12]  and  (3],  Wasow  and  Lehman  made  a study  of  the  asymptotic  behavior 
of  harmonic  functions  near  reentrant  corners.  Their  studies  were  quite  neneral, 
covering  curved  boundaries  and  a wide  variety  of  conditions.  The  series  they  got 
were  only  asymptotic,  but  series  like  (2.7)  were  typical  (except  that  (2.7) 
converges) . Indeed,  we  are  lucky  with  our  particular  problem,  in  that  our 

u*s^(x,y)  is  harmonic  in  the  neighborhood  of  the  corner.  More  generally,  u^3 * S\x,y) 
contributes  additional  complications,  such  as  terms  involving  logarithms. 


In  view  of  this,  one  wonders  why  the  authors  of  (1]  managed  to  get  such  good 
results  near  the  reentrant  corner.  This  came  about  as  follows.  In  order  to  be 
able  to  check  if  their  procedure  was  giving  the  right  answers,  they  took  a problem 
in  which  the  answers  were  known.  They  chose  u(x,y)  a function  that  was  well 
behaved  over  a much  larger  region  than  that  shown  in  Fig.  1.  From  it,  they  read 
values  around  the  boundary,  and  proceeded  to  solve,  getting  back  u(x,y)  of 
course.  Since  they  started  with  a function  that  was  well  behaved  over  a large 
region,  including  the  reentrant  corner,  they  insured  that  = 0 

in  (2.7).  So  of  course  they  had  well  behaved  fourth  derivatives,  and  (1.2)  was 
an  excellent  approximation  to  (1.1),  and  their  answers  agreed  closely  with  the 
true  values.  Very  cumforting  for  them,  but  very  misleading  for  the  reader.  Had 
they  used  the  boundary  conditions  (2.3),  their  answers  would  have  been  very  poor 
near  C.  On  the  other  hand,  they  probably  did  not  have  a way  to  get  the  correct 
answers  for  the  boundary  conditions  (2.3),  and  so  would  not  have  known  if  they 
had  good  answers  or  not. 

3.  The  solution  inside  a rectangle.  This  brings  us  to  the  crucial  ques- 

tion. How  does  one  get  correct  answers  with  boundary  conditions  like  (2.3)7 
First  we  have  to  have  a technique  for  carrying  out  a solution  inside  a rectangle, 
which  we  now  explain.  Given  a rectangle  with  smooth  boundary  conditions  prescribed 
around  its  perimeter,  how  does  one  determine  a u(x,y)  which  is  harmonic  in  the 

interior  and  takes  the  prescribed  values  on  the  perimeter? 


5 


Figure  3 


Consider 

* 

(3.1)  u (x,y)  = u(x,y)  - A - Bx  - Cy  - Dxy  . 

This  is  harmonic,  and  has  smooth  boundary  conditions.  It  is  easy  to  choose  A, 

B,  C,  and  D so  that  u (x,y)  takes  the  value  zero  at  each  corner  of  the 

rectangle.  Along  the  top  let  u (x,y)  = f (x) , for  a £ x <_  b.  Our  choice  of 

A,  B,  C,  and  D assures  that  f(a)  = f (b)  = 0.  Also,  as  we  were  assuming  smooth 
boundary  conditions,  let  us  say  that  that  assures  that  f (x)  has  almost  every- 
where a second  derivative  for  a 1.  x ;S  b which  is  of  bounded  variation.  So 
f (x)  satisfies  the  conditions  of  Lemma  2.1.  We  get  its  Fourier  expansior,  (2.1), 
with  the  D defined  by  (2.2).  Consider 


00  sinh 

ut(X'y)  = l Dm  : 


xm(y  - c) 

b - a 

, irmh 

lh  r 

b - a 


n(x  - a) 
b - a 


where  h is  the  height  of  the  rectangle  and  c is  the  value  of  y at  the 
bottom  of  the  rectangle.  Clearly  ut(x,y)  is  harmonic.  It  is  zero  along  the 

left  side  of  the  rectangle  (x  = a) , it  is  zero  along  the  right  side  of  the 
rectangle  (x  = b) , it  is  zero  along  the  bottom  of  the  rectangle  (y  * c) , and 
it  equals  f(x)  along  the  top  of  the  rectangle;  that  is,  on  the  top  it  agrees 
with  u (x,y) . We  carry  out  an  analogous  construction  for  each  of  the  other 
three  sides  of  the  rectangle,  and  add  together  the  resulting  four  series.  Since 
the  sum  agrees  with  u (x,y)  on  the  entire  perimeter,  it  has  to  be  equal  to 
u*(x,y)  throughout  the  rectangle.  Then  we  determine  u(x,y)  from  (3.1). 

Armed  with  this  technique,  let  us  return  to  the  problem  of  Fig.  1.  Values 
of  u(x,y)  have  been  prescribed  around  the  boundary  (for  example,  see  (2.3)). 

We  guess  values  along  CF.  With  the  boundary  conditions,  this  gives  values 
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around  the  perimeter  of  the  rectangle  ABFG.  As  described  just  above,  we  get  a 
harmonic  function  inside  this  rectangle  which  takes  the  prescribed  boundary  condi- 
tions. It  gives  us  values  along  HC.  With  these  and  the  boundary  conditions, 
we  have  values  around  the  perimeter  of  the  rectangle  HDEG.  From  these,  we  get 
values  in  the  interior,  including  along  the  line  CF.  This  will  be  an  improvement 
over  our  first  guess. 

We  repeat  the  process.  In  an  actual  calculation,  with  the  conditions  (2.3), 
it  took  about  fifteen  iterations  for  convergence.  However,  because  the  Fast 
Fourier  Transform  gets  the  D very  quickly,  the  calculation  to  convergence  did 
not  take  very  long.  However,  it  did  not  converge  to  the  u(x,y)  we  were  seeking. 
Recall  that  in  Lemma  2.1,  it  was  required  that  f (x)  have  a second  derivative  of 
bounded  variation.  But  the  u(x,y)  defined  by  conditions  (2.3)  has  an  infinite 
first  derivative  as  one  approaches  C along  CF. 

This  seems  too  bad.  However,  the  procedure  we  just  described  is  not  without 
value.  In  tact,  it  will  be  the  one  we  will  use  in  the  end,  but  with  a slight 
modification.  Our  difficulty  (refer  to  Fig.  1)  is  that,  along  the  lines  BF  and 
HD,  the  function  u(x,y)  that  we  are  trying  to  determine  does  not  have  almost 
everywhere  a second  derivative  of  bounded  variation.  If  we  should  try  this 
procedure  on  a u(x,y)  which  does  have  almost  everywhere  a second  derivative  of 
bounded  variation  along  the  lines  BF  and  HD,  we  would  succeed  admirably  in 
determining  that  u(x,y),  and  in  terms  of  rapidly  converging  Fourier  series.  All 
we  need  for  u(x,y)  is  to  know  its  values  around  the  boundary,  and  to  be  assured 
that  it  is  sufficiently  smooth  along  the  lines  BF  and  HD. 

4.  A slight  modification.  Recall  that  the  u(x,y)  we  are  seeking  to 
determine  equals  (2.4)  plus  {2.1)  inside  and  on  the  three-quarters  circle,  and 
that  (2.4)  is  very  smooth  along  the  lines  BF  and  HD.  Because  of  this,  we 
will  show  that 


(4.1) 


u(x,y) 


4 


has  a second  derivative  of  bounded  variation  along  both  the  lines  BF  and  HD. 
Along  BC  and  CD,  (4.1)  equals  the  right  side  of  (2.3),  which  is  very  smooth. 
Inside  and  on  the  three-quarters  circle,  (4.1)  equals  (2.4)  plus  the  remainder 
of  the  series  (2.7),  which  is  smooth  enough.  And  from  the  three-quarters  circle 
out  to  F or  H,  (4.1)  is  the  sum  of  three  harmonic  functions  out  to  a straight 
line  border  along  which  their  boundary  values  are  infinitely  differentiable;  by 
Lemma  2.2,  all  derivatives  exist  continuously  out  to  the  border. 

If  we  could  somehow  determine  the  values  of  D^  and  D^,  we  could  determine 

(4.1)  by  the  procedure  of  the  previous  section.  We  certainly  can  determine  the 
values  of  (4.1)  around  the  boundary;  we  had  had  to  choose  a value  of  A before 
the  values  of  D^  and  D^  could  be  defined  (in  fact,  we  had  chosen  A = 1 for 

our  calculation),  and  the  values  of  u(x,y)  are  given  by  (2.3).  Also,  as  we 
have  just  carefully  ascertained,  (4.1)  has  a second  derivative  of  bounded 
variation  along  the  lines  BF  and  HD.  Being  given  the  values  of  D^  and  D^. 
we  could  then  calculate  u(x,y)  from  (4.1). 
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So  we  are  faced  with  the  problem  of  determining  Dj  and  Dj. 

We  remind  the  reader  that  a computer  operates  linearly.  To  calculate  (4.1) 
by  the  procedure  of  the  previous  section,  we  would  get  the  same  numerical  answers 
by  either  of  the  two  following  procedures. 

(1)  Apply  the  procedure  to  the  total  function  (4.1). 

(2)  Apply  the  procedure  first  to  u(x.y),  getting  some  Fourier  expansions 
S1,  then  apply  the  procedure  to 


getting  some  Fourier  expansions 


then  apply  the  procedure  to 


getting  some  Fourier  expansions 
expansions  into 


and  finally  combine  the  various  Fourier 


Although  S1  will  be  a poor  representation  of  u(x,y) , as  we  observed  in 

the  previous  section,  and  S11  and  S111  will  be  poor  representations  of  (4.2) 
and  (4.3',  for  similar  reasons,  the  combination  (4.4)  will  be  an  excellent  repre- 
sentation of  (4.1),  since  the  linearity  of  the  computer  assures  that  it  comprises 
the  same  numbers  that  one  would  get  by  applying  the  procedure  of  the  previous 
section  to  the  entirety  of  (4.1). 


. I II  „ „HI 

With  no  more  ado,  let  us  proceed  to  determine  S , S , and  S 


Consid- 


ering 


D2  as  two  (as  yet)  unknown  parameters,  we  can  take  (4.4)  as 


representing  (4.1).  Subtracting  (2.4)  from  (4.4),  we  will  have  a representation 
of 

2 4 

(4.5)  u(x,y)  - u<s)(x,y)  - D^}  sin  ~ - D2(a)  sin  ^ * 

That  is,  using  (4.4)  minus  (2.4),  we  can  actually  calculate  values  of  (4.5)  at 
any  points  of  the  L-shaped  region  of  Fig.  2,  except  that  the  values  will  come 
out  as  linear  combinations  of  and  D2- 

Observe  that  (4.5)  is  zero  on  both  the  lines  BC  and  CD,  since  (2.6)  was. 
So,  by  the  same  method  that  we  used  to  get  the  expansion  (2.7)  for  (2.6)  inside 
and  on  the  three-quarters  circle,  we  can  get  an  expansion  like  (2.7)  for  (4.5). 
Obviously,  this  expansion  has  to  consist  of  (2.7)  with  the  first  two  terms 
deleted.  So,  when  we  take  m = 1 and  2 in  (2.2)  to  get  and  D^,  we  must 
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get  the  value  zero.  But,  as  the  values  of  f (x)  in  (2.2)  are  taken  from  (4.4) 
minus  (2.4),  the  numerical  quadratures  to  determine  (2.2)  must  yield  linear 
combinations  of  and  D., . Putting  these  linear  combinations  equal  to  zero 

for  m = 1 and  m = 2 gives  us  two  simultaneous  linear  equations  for  and 

• We  solve  these.  Putting  the  solutions  into  (4.4)  gives  Fourier  expansions 

for  (4.1).  But  now  we  knew  and  , and  so  can  calculate  u(x,y)  from  (4.1). 
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